Análisis Exploratorio de Suelos (Honoria - Tournavista)
Alonso Otiniano Zavala & Joseps Andrade Choque - (Dr. Jimmy Rosales)
2022-12-04
knitr::opts_chunk$set(echo = TRUE)
Sys.setenv('MAPBOX_TOKEN' = 'pk.eyJ1IjoiYWxvbnNvMjUiLCJhIjoiY2tveGJseXJpMGNmcDJ3cDhicmZwYmY3MiJ9.SbThU_R8YGE1Zll-nNrZKA')
library(readxl) # lectura archivos excel.
library(tidyverse);library(corrplot) # librearias de tratamiento de data## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## corrplot 0.92 loadedlibrary(sf); library(raster); library(rgdal) # librerias espaciales basicas## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
## Loading required package: sp
##
## Attaching package: 'raster'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/Alonso/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Alonso/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(mapview); library(mapedit); library(mapdeck) # librerias espaciales visualizacion especifica## Warning: package 'mapdeck' was built under R version 4.2.2##
## Attaching package: 'mapdeck'
##
## The following object is masked from 'package:tibble':
##
## add_columnlibrary(leaflet); library(leaflet.extras)
library(htmltools)
#ms <- mapdeck_style("satellite")
mapviewOptions(basemaps = "Esri.WorldImagery")Publicación Internacional:
Certificado de Aceptacion:
Study to Determine Levels of Cadmiumin Cocoa Crops Applied to Inland Areas of Peru
Paper:
Acreditation
Link de Descarga
El link es el siguiente: Agronomy MDPI.
Video de Presentación
El link es el siguiente: ECITEC UNI 2020
Análisis de la Información Recolectada:
Lectura de Data
Cacao <- read_xlsx(path = "data_IIC1/BD_Cd.xlsx", col_names = TRUE)
str(Cacao)## tibble [24 × 7] (S3: tbl_df/tbl/data.frame)
## $ Xeast : num [1:24] 518642 518644 518690 518701 518706 ...
## $ Ynorth : num [1:24] 9032291 9032343 9032334 9032232 9032284 ...
## $ Cdsoil_mgkg : num [1:24] 1.75 1.6 1.9 1.7 1.25 1.75 1.5 1.4 1.4 1.35 ...
## $ Pbsoil_mgkg : num [1:24] 5.41 5.03 7.24 5.81 5.91 5.62 6.57 5.26 5.01 7.02 ...
## $ Znsoil_mgkg : num [1:24] 70.5 99.4 77.1 97.1 74.7 ...
## $ pHapprox : num [1:24] 7.06 6.55 6.28 8 6.7 6.4 6.53 5.79 5.41 6.39 ...
## $ Ceapprox_uS_cm: num [1:24] 39.9 10.3 5 126.7 8 ...head(Cacao, n=5)tail(Cacao, n=5)summary(Cacao)## Xeast Ynorth Cdsoil_mgkg Pbsoil_mgkg
## Min. :518642 Min. :9031868 Min. :1.100 Min. :4.740
## 1st Qu.:519033 1st Qu.:9031913 1st Qu.:1.350 1st Qu.:5.247
## Median :519183 Median :9032036 Median :1.575 Median :5.850
## Mean :519108 Mean :9032052 Mean :1.634 Mean :5.929
## 3rd Qu.:519221 3rd Qu.:9032134 3rd Qu.:1.863 3rd Qu.:6.525
## Max. :519556 Max. :9032343 Max. :2.550 Max. :7.240
## Znsoil_mgkg pHapprox Ceapprox_uS_cm
## Min. :11.80 Min. :5.410 Min. : 4.50
## 1st Qu.:44.79 1st Qu.:6.160 1st Qu.: 9.00
## Median :58.70 Median :6.375 Median : 13.50
## Mean :58.66 Mean :6.494 Mean : 26.11
## 3rd Qu.:71.89 3rd Qu.:6.588 3rd Qu.: 40.02
## Max. :99.40 Max. :8.000 Max. :126.69# Verificar nulos
Vacios <- sapply(Cacao, function(x) sum(is.na(x)))
Vacios## Xeast Ynorth Cdsoil_mgkg Pbsoil_mgkg Znsoil_mgkg
## 0 0 0 0 0
## pHapprox Ceapprox_uS_cm
## 0 0#Con mapa:
library(Amelia)## Loading required package: Rcpp## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.8.0, built: 2021-05-26)
## ## Copyright (C) 2005-2022 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##missmap(Cacao,col=c("black","grey"),legend = FALSE)## Warning: Unknown or uninitialised column: `arguments`.
## Unknown or uninitialised column: `arguments`.## Warning: Unknown or uninitialised column: `imputations`.
# Verificar alguna condicion y eliminar valores
# id <- which(is.na(Cacao$pHapprox)) # colocar cualquier condicion
# Cacao <- Cacao[-id, ]
id <- which(Cacao$Cdsoil_mgkg>= 1.4 & Cacao$Cdsoil_mgkg <= mean(Cacao$Cdsoil_mgkg, na.rm=TRUE))
Cacaonew <- Cacao[id, ]
Cacaonew %>% dplyr::select(Pbsoil_mgkg) %>% summarise(media_Pb = mean(Pbsoil_mgkg))Exploracion Basica
Univariante
par(mfrow=c(2,2))
hist(Cacao$Cdsoil_mgkg)
boxplot(Cacao$Cdsoil_mgkg)
plot(density(Cacao$Cdsoil_mgkg))
plot(log(Cacao$Cdsoil_mgkg))
Revisión del Modelo de Elevación Digital:
DEM <- raster("data_IIC1/DEM_DRONE/Zona_Oeste.tif")
str(DEM)## Formal class 'RasterLayer' [package "raster"] with 12 slots
## ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots
## .. .. ..@ name : chr "D:\\3.0Adicionales\\FIGMM_SMAUML\\ParteII\\SlidesR3\\data_IIC1\\DEM_DRONE\\Zona_Oeste.tif"
## .. .. ..@ datanotation: chr "INT1U"
## .. .. ..@ byteorder : chr "little"
## .. .. ..@ nodatavalue : num -Inf
## .. .. ..@ NAchanged : logi FALSE
## .. .. ..@ nbands : int 4
## .. .. ..@ bandorder : chr "BIL"
## .. .. ..@ offset : int 0
## .. .. ..@ toptobottom : logi TRUE
## .. .. ..@ blockrows : int [1:4] 256 256 256 256
## .. .. ..@ blockcols : int [1:4] 256 256 256 256
## .. .. ..@ driver : chr "gdal"
## .. .. ..@ open : logi FALSE
## ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots
## .. .. ..@ values : logi(0)
## .. .. ..@ offset : num 0
## .. .. ..@ gain : num 1
## .. .. ..@ inmemory : logi FALSE
## .. .. ..@ fromdisk : logi TRUE
## .. .. ..@ isfactor : logi FALSE
## .. .. ..@ attributes: list()
## .. .. ..@ haveminmax: logi FALSE
## .. .. ..@ min : logi NA
## .. .. ..@ max : logi NA
## .. .. ..@ band : int 1
## .. .. ..@ unit : chr ""
## .. .. ..@ names : chr ""
## ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots
## .. .. ..@ type : chr(0)
## .. .. ..@ values : logi(0)
## .. .. ..@ color : logi(0)
## .. .. ..@ names : logi(0)
## .. .. ..@ colortable: logi(0)
## ..@ title : chr(0)
## ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots
## .. .. ..@ xmin: num 518578
## .. .. ..@ xmax: num 519198
## .. .. ..@ ymin: num 9031751
## .. .. ..@ ymax: num 9032476
## ..@ rotated : logi FALSE
## ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
## .. .. ..@ geotrans: num(0)
## .. .. ..@ transfun:function ()
## ..@ ncols : int 24794
## ..@ nrows : int 28989
## ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs"
## .. .. ..$ comment: chr "PROJCRS[\"WGS 84 / UTM zone 18S\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\"| __truncated__
## ..@ history : list()
## ..@ z : list()# DEM
# slope <- terrain(DEM, opt ="slope")
# aspect <- terrain(DEM, opt="aspect")
# DEM.hill <- hillShade(slope, aspect, angle = 40, direction = 270)
# esta es importante pero demora DEM <- crop(DEM, extent(518600, 519200, 9031800, 9032400))
#plot(DEM, main = "Drone - Digital Elevation Model (DEM)", col = topo.colors(20),
# zlim=c(0, 250))
image(DEM, main = "Drone - Digital Elevation Model (DEM)", col = topo.colors(20),
zlim=c(0, 250))
Revisión del DEM con factores:
Factores_cacao <- st_read(dsn = "data_IIC1/Cacao_Factors/Age_Planta.shp")## Reading layer `Age_Planta' from data source
## `D:\3.0Adicionales\FIGMM_SMAUML\ParteII\SlidesR3\data_IIC1\Cacao_Factors\Age_Planta.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 518541 ymin: 9031841 xmax: 519584 ymax: 9032609
## Projected CRS: WGS 84 / UTM zone 18Sstr(Factores_cacao)## Classes 'sf' and 'data.frame': 6 obs. of 6 variables:
## $ Id : num 0 1 2 3 4 5
## $ Age : chr "E1" "E5" "E6" "E4" ...
## $ Shape_Leng: num 1085 1103 806 1199 1108 ...
## $ Shape_Area: num 69566 73596 38867 66861 51878 ...
## $ Uso_Tierra: chr "Pasto 10" "Descanso 3" "Maiz 10" "Maiz 10" ...
## $ geometry :sfc_POLYGON of length 6; first list element: List of 1
## ..$ : num [1:7, 1:2] 519345 519313 519527 519584 519525 ...
## ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
## ..- attr(*, "names")= chr [1:5] "Id" "Age" "Shape_Leng" "Shape_Area" ...summary(Factores_cacao)## Id Age Shape_Leng Shape_Area
## Min. :0.00 Length:6 Min. : 805.5 Min. :38867
## 1st Qu.:1.25 Class :character 1st Qu.:1084.5 1st Qu.:52966
## Median :2.50 Mode :character Median :1093.8 Median :61545
## Mean :2.50 Mean :1064.0 Mean :59500
## 3rd Qu.:3.75 3rd Qu.:1106.5 3rd Qu.:68890
## Max. :5.00 Max. :1198.6 Max. :73596
## Uso_Tierra geometry
## Length:6 POLYGON :6
## Class :character epsg:32718 :0
## Mode :character +proj=utm ...:0
##
##
## Factores_cacao <- Factores_cacao %>% mutate_if(is.character, as.factor)
summary(Factores_cacao)## Id Age Shape_Leng Shape_Area Uso_Tierra
## Min. :0.00 E1:1 Min. : 805.5 Min. :38867 Descanso 3 :1
## 1st Qu.:1.25 E2:1 1st Qu.:1084.5 1st Qu.:52966 Maiz 10 :3
## Median :2.50 E3:1 Median :1093.8 Median :61545 Maiz 7 Descanso 3:1
## Mean :2.50 E4:1 Mean :1064.0 Mean :59500 Pasto 10 :1
## 3rd Qu.:3.75 E5:1 3rd Qu.:1106.5 3rd Qu.:68890
## Max. :5.00 E6:1 Max. :1198.6 Max. :73596
## geometry
## POLYGON :6
## epsg:32718 :0
## +proj=utm ...:0
##
##
## mapview(Factores_cacao)mapview(DEM, legend =TRUE)+
mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10,
lwd = 2, color = "blue")## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 718753266
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 718753266 '## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definitionGenerando Shapefiles (Opcional):
Factores_cacao2 <- as.data.frame(Factores_cacao)
#Factores_cacao2$Uso_Tierra <- edit(Factores_cacao2$Uso_Tierra)
Factores_cacao2$Var_cocoa <- c("Creole-Aromatic","CCN51","","CCN51","Creole-Aromatic",
"CCN51")
Factores_cacao3 <- st_as_sf(Factores_cacao2)
mapview(Factores_cacao3, zcol="Var_cocoa")Honoria_elementos <- readxl::read_xlsx(path = "data_IIC1/BD_Cd.xlsx")
str(Honoria_elementos)## tibble [24 × 7] (S3: tbl_df/tbl/data.frame)
## $ Xeast : num [1:24] 518642 518644 518690 518701 518706 ...
## $ Ynorth : num [1:24] 9032291 9032343 9032334 9032232 9032284 ...
## $ Cdsoil_mgkg : num [1:24] 1.75 1.6 1.9 1.7 1.25 1.75 1.5 1.4 1.4 1.35 ...
## $ Pbsoil_mgkg : num [1:24] 5.41 5.03 7.24 5.81 5.91 5.62 6.57 5.26 5.01 7.02 ...
## $ Znsoil_mgkg : num [1:24] 70.5 99.4 77.1 97.1 74.7 ...
## $ pHapprox : num [1:24] 7.06 6.55 6.28 8 6.7 6.4 6.53 5.79 5.41 6.39 ...
## $ Ceapprox_uS_cm: num [1:24] 39.9 10.3 5 126.7 8 ...Honoria_elementos <- Honoria_elementos %>% rename(Este= Xeast, Norte=Ynorth,
Cd = Cdsoil_mgkg,
Pb = Pbsoil_mgkg,
Zn = Znsoil_mgkg,
pH = pHapprox,
Ce = Ceapprox_uS_cm)
Honoria_elementos2 <- Honoria_elementos[ ,c("Norte","Este")]
Honoria_elementos2 <- Honoria_elementos2[ ,order(c(names(Honoria_elementos2)))]
sputm <- SpatialPoints(Honoria_elementos2, proj4string=CRS("+proj=utm +zone=18 +south +datum=WGS84"))
spgeo <- spTransform(sputm, CRS("+proj=longlat +datum=WGS84"))
spgeo <- as.data.frame(spgeo)
colnames(spgeo) <- c("lng","lat")
Honoria_elementos2 <- cbind(Honoria_elementos,spgeo)
Honoria_elementos3 <- st_as_sf(Honoria_elementos2, coords = c("lng", "lat"),
remove = FALSE, crs = 4326, agr = "constant")
mapview(Honoria_elementos3, zcol = "Cd") +mapview(Factores_cacao3, zcol="Var_cocoa")+
mapview(Factores_cacao3, zcol = "Uso_Tierra")+
mapview(Factores_cacao3, zcol = "Age")Factores_cacao4 <- Factores_cacao3 %>% st_transform(crs = 4326 )
H <- Honoria_elementos3 %>% st_intersection(Factores_cacao4)
mapview(H, zcol = "Uso_Tierra", cex = "Cd")+
mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10,
lwd = 2, color = "blue")mapview(H, zcol = "Age", cex = "Cd") +
mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, lwd = 2, color = "blue")mapview(H, zcol = "Var_cocoa", cex = "Cd") +
mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, lwd = 2, color = "blue")Estandar de Calidad de Cd en Suelos:
Valor límite del Cd en suelos es: \[Cd<=1.4\]
mapview(Honoria_elementos3, zcol = "Cd", at = c(0,1.4,2.5), legend = TRUE,
col.regions = c("green", "red"))+
mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10,
lwd = 2, color = "blue")Data Table:
Factores_cacao_n <- Factores_cacao
Honoria_elementos_n <- Honoria_elementos3 %>% st_transform(crs = st_crs(Factores_cacao))
# sf_use_s2(TRUE) #cambiar a spherical geometry
intersection <- st_intersection(Honoria_elementos_n, Factores_cacao_n)
intersectionHonoria_statical <- intersection %>% st_set_geometry(NULL)Statical Summary:
min_max_mean_sd <- list(
min = ~min(.x, na.rm = TRUE),
max = ~max(.x, na.rm = TRUE),
mean = ~mean(.x, na.rm = TRUE),
sd = ~sd(.x, na.rm = TRUE)
)Honoria_statical %>%
summarise(across(c(Cd, Pb, Zn, pH, Ce), min_max_mean_sd))Honoria_statical %>% group_by(Age) %>%
summarise(across(c(Cd, Pb, Zn, pH, Ce), min_max_mean_sd))Ver vignette("colwise") para más detalles.
Boxplot Analysis:
Cd total:
ggplot(data = Honoria_statical)+
geom_boxplot(mapping = aes(x = Cd))+
coord_flip()
Cd by Factor Age and Land Use:
ggplot(data = Honoria_statical)+
geom_boxplot(mapping = aes(x = Cd))+
facet_grid(. ~ Age)+
coord_flip()
ggplot(data = Honoria_statical)+
geom_boxplot(mapping = aes(x = Cd))+
facet_grid(. ~ Uso_Tierra)+
coord_flip()
Cd by Prior Exploration Field:
ggplot(data = Honoria_statical)+
geom_boxplot(mapping = aes(x = Cd))+
facet_grid(. ~ Uso_Tierra)+
coord_flip()
Average Cadmiunm Soil
Corrplot Bivariant
correlacion <- cor(Honoria_statical[ ,c("Cd", "Pb", "Zn", "pH", "Ce")]
, method = "pearson")
corrplot(correlacion, method ="circle", diag = TRUE, title= "Correlacion de Pearson Honoria",
hclust.method = "median", addrect = 2)
corrplot(correlacion, method ="circle", diag = TRUE, title= "Correlacion de Pearson Honoria")
correlacion2 <- cor(Honoria_statical[ Honoria_statical$Uso_Tierra=="Maiz 10", c("Cd", "Pb", "Zn", "pH", "Ce")], method = "pearson")
# correlacion2 <- cor(Honoria_statical %>% filter(Uso_Tierra=="Maiz 10") %>%
# select(Cd, Pb, Zn, pH, Ce), method = "pearson")
corrplot(correlacion2, method ="number")
corrplot(correlacion2, method ="circle", order = "AOE")
corrplot(correlacion2, method = 'shade', order = 'AOE', diag = FALSE) 
corrplot(correlacion2, method = 'square', order = 'AOE', diag = FALSE, type = "lower")
correlacion3 <- cor(Honoria_statical[Honoria_statical$Uso_Tierra=="Pasto 10", c("Cd", "Pb", "Zn", "pH", "Ce")], method = "pearson")
corrplot(correlacion3, method ="circle")
corrplot.mixed(correlacion3, order = 'AOE')
Para ver toda la funcionalidad revisar : An Introduction to corrplot Package.
Opcional
Biplot Multivariant (Principal Component Analysis)
Cadmium Final Geostatistics Map:
library(gstat)class(Honoria_elementos3)## [1] "sf" "data.frame"summary(Honoria_elementos3)## Este Norte Cd Pb
## Min. :518642 Min. :9031868 Min. :1.100 Min. :4.740
## 1st Qu.:519033 1st Qu.:9031913 1st Qu.:1.350 1st Qu.:5.247
## Median :519183 Median :9032036 Median :1.575 Median :5.850
## Mean :519108 Mean :9032052 Mean :1.634 Mean :5.929
## 3rd Qu.:519221 3rd Qu.:9032134 3rd Qu.:1.863 3rd Qu.:6.525
## Max. :519556 Max. :9032343 Max. :2.550 Max. :7.240
## Zn pH Ce lng
## Min. :11.80 Min. :5.410 Min. : 4.50 Min. :-74.83
## 1st Qu.:44.79 1st Qu.:6.160 1st Qu.: 9.00 1st Qu.:-74.83
## Median :58.70 Median :6.375 Median : 13.50 Median :-74.83
## Mean :58.66 Mean :6.494 Mean : 26.11 Mean :-74.83
## 3rd Qu.:71.89 3rd Qu.:6.588 3rd Qu.: 40.02 3rd Qu.:-74.83
## Max. :99.40 Max. :8.000 Max. :126.69 Max. :-74.82
## lat geometry
## Min. :-8.758 POINT :24
## 1st Qu.:-8.758 epsg:4326 : 0
## Median :-8.757 +proj=long...: 0
## Mean :-8.757
## 3rd Qu.:-8.756
## Max. :-8.754Bubbgle Map:
Honoria_elementos4 <- Honoria_elementos3 %>% st_set_geometry(NULL)
coordinates(Honoria_elementos4) <- ~Este+Norte
class(Honoria_elementos4)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"summary(Honoria_elementos4)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## Este 518642 519556
## Norte 9031868 9032343
## Is projected: NA
## proj4string : [NA]
## Number of points: 24
## Data attributes:
## Cd Pb Zn pH
## Min. :1.100 Min. :4.740 Min. :11.80 Min. :5.410
## 1st Qu.:1.350 1st Qu.:5.247 1st Qu.:44.79 1st Qu.:6.160
## Median :1.575 Median :5.850 Median :58.70 Median :6.375
## Mean :1.634 Mean :5.929 Mean :58.66 Mean :6.494
## 3rd Qu.:1.863 3rd Qu.:6.525 3rd Qu.:71.89 3rd Qu.:6.588
## Max. :2.550 Max. :7.240 Max. :99.40 Max. :8.000
## Ce lng lat
## Min. : 4.50 Min. :-74.83 Min. :-8.758
## 1st Qu.: 9.00 1st Qu.:-74.83 1st Qu.:-8.758
## Median : 13.50 Median :-74.83 Median :-8.757
## Mean : 26.11 Mean :-74.83 Mean :-8.757
## 3rd Qu.: 40.02 3rd Qu.:-74.83 3rd Qu.:-8.756
## Max. :126.69 Max. :-74.82 Max. :-8.754bubble(Honoria_elementos4, "Cd", col = c("green","green"), main = "Cadmium Concentrations (ppm)")
Grilla generada por Meuse Data:
data(meuse.grid)
ggplot(data = meuse.grid) + geom_point(aes(x, y))
# Hacemos una grilla regular con `SpatialPolygonsDataFrame`
proj4string(Honoria_elementos4) <- CRS("+proj=utm +zone=18 +south +datum=WGS84") #definimos proyecciónHonoria_elementos4 <- spTransform(Honoria_elementos4, CRS("+proj=longlat +datum=WGS84"))
grd <- makegrid(Honoria_elementos4, n = 1000) #generamos grilla
colnames(grd) <- c("x", "y") #asignamos nombres
plot(grd$x,grd$y) #visualizamos
# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
coords = grd,
proj4string = CRS(proj4string(Honoria_elementos4))
)
# Then, visualize your clipped grid which can be used for kriging
ggplot(as.data.frame(coordinates(grd_pts))) +
geom_point(aes(x, y))
gridded(grd_pts) <- TRUE
grd_pts <- as(grd_pts, "SpatialPixels")Limite2 <- readOGR(dsn="data_IIC1/Cacao_shp//Limite_Cacao..shp")## OGR data source with driver: ESRI Shapefile
## Source: "D:\3.0Adicionales\FIGMM_SMAUML\ParteII\SlidesR3\data_IIC1\Cacao_shp\Limite_Cacao..shp", layer: "Limite_Cacao."
## with 1 features
## It has 1 fields
## Integer64 fields read as strings: Idgrd <- makegrid(Limite2, n = 1000) #generamos grilla
colnames(grd) <- c("x", "y") #asignamos nombres
plot(grd$x,grd$y) #visualizamos
# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
coords = grd,
proj4string = CRS(proj4string(Limite2))
)
# subset all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[Limite2, ]
grd_pts_in2 <- grd_pts_in
proj4string(grd_pts_in2) <- CRS("+proj=utm +zone=18 +south +datum=WGS84")
gridded(grd_pts_in2) <- TRUE
grd_pts_in2 <- as(grd_pts_in2, "SpatialPixels")cd.idw = idw(Cd~1, Honoria_elementos4, grd_pts_in)## [inverse distance weighted interpolation]spplot(cd.idw["var1.pred"], main = " Cadmium inverse distance weighted interpolations")
cd.idw = idw(Cd~1, Honoria_elementos4, grd_pts_in2)## [inverse distance weighted interpolation]spplot(cd.idw["var1.pred"], main = " Cadmium inverse distance weighted interpolations")
m <- vgm(.59, "Sph", 874, .04) #setear bien
x <- krige(log(Cd)~1, Honoria_elementos4, grd_pts_in2, model = m)## [using ordinary kriging]spplot(x["var1.pred"], main = "ordinary kriging predictions")
spplot(x["var1.var"], main = "ordinary kriging variance")